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Gntzwiller wavefunction is a physically well motivated trial wavefunction for describing 
correlated electron systems. In this work, a new approximation is introduced to facilitate 
evaluation of the expectation value of any operator within the Gntzwiller wavefunction for¬ 
malism. The basic idea is to make use of a specially designed average over Gntzwiller 
wavefunction coefficients expanded in the many-body Fock space to approximate the ratio 
of expectation values between a Gntzwiller wavefunction and its underlying noninteracting 
wavefunction. To check with the standard Gntzwiller approximation (GA), we test its per¬ 
formance on single band systems and find quite interesting properties. On finite systems, we 
noticed that it gives superior performance than GA, while on infinite systems it asymptoti¬ 
cally approaches GA. Analytic analysis together with numerical tests are provided to support 
this claimed asymptotical behavior. At the end, possible improvements on the approximation 
and its generalization towards multiband systems are illustrated and discussed. 
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INTRODUCTION 


Correlation effects play an important role in electronic movements and physical properties of 
real materials[l|. Strong electron-electron interaction is believed jto be key to fully understand 


many interesting phenomena, including high Tc superconductivity!^, heavy fermion behaviors 
abnormal transport and optical propertiesjj], and more on strongly correlated materials. A versa¬ 
tile and convenient way to describe these systems has been posed as a big challenge to the solid state 
community in the past years. Part of the reasons are due to the strong Coulomb interaction pre¬ 
venting a controlled perturbative treatment which has been very successfully developed for weakly 
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interacting systems. The conventional mean field treatment is most versatile in gaining insights to 
correlated electron systems on their possible new physics, competing phases and dynamical behav¬ 
iors. But its conclusions are always in question whether the ignored residual electron correlation 
effects are still significant enough to overshadow the presumed mean-field behaviorsjsjjg]. To go 
beyond the mean field approximation, analytic tools were proposed based on infinite summation 
of Feynmann diagrams of chosen types, e.g., the Random Phase Approximation (RPA) T]], the 
spin fluctuation theory J ], and compact diagrammatic equations like the Fluctuation Exchange 
approximation (FLEX) 9| or the parquet formalism [l^ [ll |. Although they offer alternative ways 
for people to use, their results might be biased by starting with a subjective choice of a subgroup 
of diagrams which are usually numerically convenient to deal with. On the side of computational 


physics, a number of tools exist including Exact Diagonalization (ED)[l^, Quantum A 


simulation (QMC) [13|] , Dynamic Mean Eie’ 


d Theory (DMET) and its cluster extensions 


and renormalization group type methods 17l | 


]QQ- 


onte Carlo 


These tools are computationally demanding 


and can suffer from serious finite size effects or other numerical complications. 

On the other hand, the variational approachhas always been a very important category of 
methods in addressing a wide range of physical problems, thanks to its capability of conveniently 
incorporating clear physical insights into a well-designed trial wavefunction. A typical example 
is the ab initio local density approximation (LDA) which has become an indispensable tool in 


modern scientific and material research [211] 
approach can play an important role and 
studying various correlation effects [23|] [24] 


. In strongly correlated systems, the variational 
ras been shown to be very effective and enlightening in 


25(1 . It avoids looking for a small quantity to validate a 


perturbative expansion of the problem, but instead, it focuses directly on most prominent physics 
out of correlation effects. A simple-minded way to take care of the local onsite correlation is to 
form a trial wavefunction with a local projection operator to control double occupancy. This brings 
up the famous Gutzwiller trial wavefunction to address local correlation effects 25(|. Unfortunately, 
a direct evaluation of operator expectation values is not practical for a many-body wavefunction. 
How to form a new way to efficiently implement the Gutzwiller wavefunction formalism will be the 
focus of the current paper. 

The typical model Hamiltonian to study correlation effects is the so-called single band Hubbard 


model 


261 ] on a two-dimensional square lattice, expressed as 


H = -t ^ 

dJ>w i 


( 1 ) 


with i,j denoting site indices and a spins. Two energy scales, the hopping amplitude, t, between 
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nearest neighboring sites and the onsite Hubbard interaction, [/, exist and compete with each 
other. The correlation effect becomes dominant in case onsite energy U is big as compared against 
hopping t. The ground state wavefunction is presumably well described by the Gutzwiller trial 
wavefunction (GWF) defined as 

with the local Gutzwiller projector defined as 

= (3) 


Here h-j-hj, is defined as the local double occupancy operator and g is the unknown Gutzwiller 


variational parameter, a positive number which controls the weight of e 
containing double occupancy in the noninteracting wavefunction ITq 


ectronic configurations 


25!]. Although Pg is 


defined purely locally, it might still capture some nonlocal physics through the 2nd order virtual 
hopping process in the strong correlation limit 


n 


281 ] [29[ . Thus Eq. [2] is a quite reasonable trial 


wavefunction to capture the essential physics of a correlated electron system. 

However, introducing a physically sensible trial wavefunction is only the first step. Brute force 
evaluation of physical observables in a many-body system can be a real pain due to numerical 


difficulties and computer hardware capacity. 


was introduced to facilitate the evaluation 


rus, t.^ 


re so-called Gutzwiller approximation (GA) 


25 ! aqlsii d. There, 


the expectation value for an 


interacting system is assumed to be proportional to that of the noninteracting wavefunction through 
a site-decomposeable renormalization factor. Taking the hopping operator as an example, GA gives 




( 4 ) 


for two distinct sites. Here, a, fi denote orbital indices and a denotes spin. A succinct analytic 
expression for Zi^a is available for the operator (or Ci^a,a) 32 j. Actually, we might come up 


with an educated guess that, given any local operator 6/ with I the set of local spin-orbital states 
involved in d, the corresponding zq is given as 


1 


Y^n(a,o-)G/ ”(*a,(T) (1 ”(*a,(T)) 


'^Vprpr |(r'| 6/ |r) 


( 5 ) 


with 


A 


(a,a) 


= 


n 


-— if h(a ct) is part of dj 

1 — nV N ^ ^ 

{a,a) 

1 o.w. 


( 6 ) 
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Here iT'^aa) ~ ('^^ol ^(a,o-) l^o) and pr denotes the occupational probability at a local Fock state 
r. The summation in Eq. [5] traverses all local Fock states mutually related by the local operator 
d/. One can go through the tedious yet well-established procedure to verify this relations for any 
operator 33|]. The introduction of Zi^a as a renormalization factor is physically very sensible. Strong 


correlation effects among electrons necessarily leave an impact on their motions, thus modifying 
their dynamical and transport behaviors as compared against a Fermi liquid system described by 

l^o). 

There have been roughly three different ways to reach GA. Metzner et al approached GA with 


mathematica 

wavefunction 


ri gor 


by applying Feymann diagrammatic expansion techinque to the Gutzwiller 


34i | 35l | . The approach effectively cuts intersite communications in the infinite spatial 


dimension limit and leaves only intrasite correlations among local orbits. This brings up the 
main site-wise decomposibility feature of GA readily seen from the site-dependent renormalization 
factors illustrated in Eq. HI This approach has been further carried out by Bunemann et. al. and 
arizio et. al. to multiband systems and with a more generalized Gutzwiller projector 


37 


3^[36| 


38l |. The second way of writing down GA is through physical intuition with a hand-waving 


argument. Gutzwiller originally formulated GA based on some assumptions 


25l | whose physics were 


rather obscure thus preventing it from being generalized to multiband systems. Ogawa’s counting 
argument formally pointed out that the physics underlying Gutzwiller’s GA is to assume that the 
expectation value of a product of number operators equals the product of expectation value of 


each number operator in the series 


3ll |. The physics of GA is most transparent in Bunemann’s 


version of counting argument [3^. We regard it as a third approach towards GA due to its balance 
between clarity in physics and completeness in formulation. Bunemann made it clear that the 
projection of a noninteracting wavefunction onto a specific electronic configuration depends only 
on the number of electrons on each local orbital on the lattice, but not on how these electrons are 
distributed on the lattice. This implies two things. Eirst, the combinatoric trick accompanying the 
assumption can be used to evaluate any inner product, which greatly simplifies the whole formalism 
and introduces an additional convenience of taking the thermodynamic limit to further simplify 
the expressions. Second, he pointed out that GA actually implies more than just one relation, but 
instead, infinitely many for different operators whose expectation values are to be evaluated. This 
leaves us a taste on how crude GA might work in reality. 

Despite many achievements made by GA on qualitatively addressing the Mott physics in cor¬ 
related electron systems, quantitatively however, it might introduce artifacts. For instance, GA 
predicts existence of a Brinkman-Rice metal-insulator transition at a critical onsite interaction 
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comparable to the bandwidth in a single band system 


give metallic behavior unless U ^ oo 


391], while GWF formalism should always 


34l |. These artifacts are of course closely related with the 


rough assumptions made in GA. The artificial Brinkman-Rice transition seems to imply excessive 
local correlations were introduced by GA. In addition, besides the unphysical infinite spatial di¬ 
mension limit, the most questionable assumption in GA is that it calls for a presumably known 
connection between the Gutzwiller local orbital occupations and those of the noninteractiM wave- 


function, a requisite to facilitate the Feynman diagrammatic expansion to derive GA 


^[331]. This 


might cause trouble in multiband systems where charge flow among orbitals might necessarily play 
a role due to correlation effects fl¬ 
it would be useful if some of the concerns mentioned above could be addressed. However, 
within the existing GA formulation, this is not an easy task mainly because the formulation is 
mathematically prohibitively complicated to be further improved. Thus we need a new way to 
look into the renormalization factors of operators, central quantity of GA to bridge between the 
Gutzwiller and the noninteracting wavefunction. We come up with a very simple way to design 
the renormalization factors which gives the standard GA in the thermodynamic limit, and which 
also has the potential to be improved systematically. We might thus call it a fourth way to reach 
GA. 


METHOD 

Let’s introduce some notations first to facilitate the up-coming discussion. Given a number of 
local spin-orbital states, the local Fock space is set up and denoted as F. Specifically, in the single 
band case, we denote the spin-orbital states as 1 = (a,t), 2 = with a the orbital index, and 

note F E {0, (1), (2), (I, 2)}. Here 0 denotes an empty state and (1) denotes a singly occupied state 
with spin-orbital state I taken. Similar interpretation applies to (2) and (1,2). Let’s denote jF] to 
represent the number of electrons in a Fock state F, and Sz (F) the total spin z component of F. 
Specifically, we have the following 


|0|=O,S. (0) = O; 

|(1)| = 1.S.(1) = I; 

|(2)| = 1,«.(2) = -I; 
|( 1 . 2 )| = 2 ,«.( 1 , 2 ) =0 
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Suppose the system has a total of N sites, electrons and net spins. Each electronic con¬ 
figuration is described in an occupation representation as {E} = {ri,r2 ,... jEat} with subscripts 
denoting site indices. These occupation representations form a complete Fock space on the lattice, 
satisfying 

Y^\r){T\=i ( 7 ) 

{r} 

The count of each local Fock state in E, nr (E) , is defined to be 

nr (E) = ^ m^r (E) with n^^r (E) = <5r,ri (8) 

i 

or simply nr if no confusion arouses in its interpretation. All possible values for nr form a set 
denoted as {nr} . 

Here is something quite general and useful for us to know. Consider a generic trial wavefunction 
|T) and note the closure relation of Fq. [3 the inner product involving any operator O is written 
as 


(^|0 |^)= (^|r) (E|0|E') (E'|^> 


(9) 


{r},{r'} 

This can be reexpressed in terms of a noninteracting wavefunction |To) on which IT) is set up as 


with 


(T|d|T)= ^ C(E,E') (To|E)(E|d|E')(E'|To) 


( 10 ) 


r(r WlilW 

^ ^ (To|E) (E'lTo) 


( 11 ) 


which holds as long as the denominators do not vanish, a fairly mild constraint to be met in most 
cases. If we manage to replace C (E, E') with some constant C, an idea borrowed from the first 
mean value theorem of integrals, Fq. [TU]is now 


(T|0 |T)~C ^ (To|E)(E|0|E')(E'|To) 

= C(To|0|To) 


( 12 ) 


In the like fashion, the expectation value of O can be related to that of the noninteracting system 
through 


/^\_ (T|d|T) ^ .T (To|0|To) _ ^ 

\ / (T|T) B (To|To) /o 


(13) 
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with 


o - B 


(14) 


Here J- and B are specifically chosen symbols to denote the constant prefactors in the numerator and 

32t |. we call the renormalization factor for operator 


the denominator respectively. Jnst like GA 


O. Thus, through a renormalization factor, the expectation value of an operator evaluated with a 
correlated Gutzwiller wavefunction is directly related to that of a noninteracting wavefunction. In 
principle, Eq. [13] is able to rigorously hold if J- and B are chosen correctly. However, this is most 
likely not the case in practice. How well Eq. [13] holds depends on the specihc method used to set 
up Z^. 

We believe the implication of the form of a renormalization factor in Eq. [13] reaches far beyond 
GA in that it gives us a clear route to determine Zq without resorting to the very involved algebras 
from Feynmann diagrams and the unphysical assumptions used in GA. One way to determine Z^, 
or equivalently T and H , is to go through C'(r,r') directly by averaging these wavefunction 
coefficients in the many-body Fock space. You will see that, in the single band case, even a simple 
dehnition towards J- and B through a wavefunction coefficient average readily gives GA as its 
limiting behavior. Here are the details for dehning T and B. 

For the Gutzwiller trial wavefunction dehned in Eq. [2] and Eq. [3] with only density operators 
entering the Gutzwiller projector, one possible way to dehne C for any operator O is 


Co {g) = 


E{r},{r'} 



5 

(o;r,r') 



fo;r,F) 



(15) 


where (r) = with To = c|c| |0) denoting the doubly occupied Fock state. The delta 

functions, 


4"^ = M - iV ] ^ ( |r| - iVe ] -5 ( (T) - Sz 


(16) 


\r /\r /\r / 

ensure conservation of total charge and spin, and expansion on all sites. The indicator function 
S ("0;r,r') is defined as 


5 0 ;r,r' = 


1 if (r|o|r) / 0 

0 o.w. 


(17) 


Thus the operator O has its effect included through S T, F'j . Although looked complicated, Eq. 

[TSlis just a carefully designed average over nonvanishing elements of the inner product 
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expanded in the Fock space set up on the whole lattice. Thus it readily guarantees the correct 
limiting behavior as g approaches unity. Why Cq can be conveniently written out in Eq. [15] is 
closely related to the current choice of the Gutzwiller projector defined on number operators only. 
An additional advantage with such a Gutzwiller projector is, this enables us to use combinatorics 
to simplify the whole expression. 

Now we are ready to define and its renormalization factor Zq. Without loss of generality, 
let’s assume O = Oijki acts on a sublattice TZ = {i,j, k, 1} defined on the affected sites, and denote 
the number of distinct sites in TZ to be Nti. We define TZ as the complementary lattice to TZ and 
denote its number of sites to be N^. A Fock space on the whole lattice can be decomposed into 
Fock spaces on the sublattices TZ and iz. The counts on local Fock states in both sublattices are 
denoted as nr and nr respectively. Then, apply Eq. [15] to both numerator and denominator and 
we have 


( 


Oijki 


) 




T 

B 



0 


(18) 


with 


^=Co (g) 


B=Ci{g) 


E'm.in {(n,«5 

'o;r,r) 


} 

E'|r),(r){'S(o;r,r') 

E’p,,cS5>.5(3) 

} 




(19) 


( 20 ) 


Here I denotes the identity operator. The primed sum in Eq. [19] denotes that {F} and {F'} are 
confined within sublattice TZ only. The double primed sum assumes hr and hr^, are counted within 
sublattice R for fixed F. The delta functions states the conservation constraints reenforced on 
TZ are 


m = ^hr - U J^hr |F| - w <5 J^hrs, (F) - Sz 


( 21 ) 


with 


Ne = Ne- Nf (F) 
Sz = Sz- sf (F) 


( 22 ) 

(23) 


where (F) and (F) denote total number of electrons and net spins on sublattice TZ occu¬ 
pied by a specific electronic configuration F. Note, a physical operator would not change electron 
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occupation and total spin z component, which implies (r) = (F') and 5^ (F) = (F'). 

takes the standard definition of the multinomial coefficient, 


^ N\ 

^ Or (^r!) 


(24) 


Eq. [18] to Eq. [23] completes the expressions for approximating the renormalization factors through 
wavefunction coefficient averaging. These expressions can be directly used to approximate the 
energy expectation value of a Hamiltonian of a small system. 

For big systems, one needs to take the thermodynamic limit to simplify the coefficients in Eq. 
[TSl A typical sum involved in Eq. [19] and Eq. [20] is 


^({nr})= (25) 

{%} 

In the thermodynamic limit, the sum can be replaced by a single term located at some unknown 
{hp} where the terms are peaked. To solve for {hp} we introduce Lagrange multipliers to relax 
the constraints entering Eq. [25] and replace the factorials with their asymptotic expressions using 
the Sterling’s formula. We finally come up with the following functional 


h ({hr} ,a,/3,^) = N In N — hp In hr + E 2 hr In^r 

r r 

+ a hr - iV^ + hr |r| - + 7 hrs, (T) - 5,^ (26) 

whose minimization gives {hp} . The Gutzwiller variational parameter, g, now carries a local Fock 
state index for the purpose of formalism consistency. Following the standard procedure to extremize 
h ({hr} , a, (3, 7 ) by taking partial derivatives with respect to {hr}, we get explicit expressions for 
hp as 

^ ^ ffpexp (/ 3 |r| +7s^ (T)) 

N Er5rexp(^|r|+7s^ (T)) 

The remaining unknown, (3 and 7 , satisfy the following two nonlinear equations, 


^ (|r| - NJn) gl exp (/3 |r| + 7 s, (T)) = 0 (28) 

r 

( 5 . (E) - ~SJn) gl exp (/3 |r| + 7 s, (T)) = 0 (29) 

r 

Note iVe and Sz defined in Eq. [2^ and Eq. [23] explicitly enter the above two equations. This implies 
that each term of the primed sums in 3F and B will have its own optimal {hp}. This complicates 
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the whole computation, but is expected to have a mild impact on its performance as Eq. [28] and 
Eq. |29| are usually well-behaved. 

Eq. [18] and its accompanying definition of Eq. [19] and Eq. [20] looks fundamentally different 
from GA in both the underlying assumptions and their mathematical formulations. Surprisingly, 
however, the current scheme can be shown to be well related to GA in the single band system. 
A brief proof is provided in Appendix A and a numerical study is given in Fig. [4] in the main 
text. It clearly reveals that GA can be regarded as being an overly simplihed approximation of 
the current scheme, but with vanishing difference at infinite lattice size. Thus, we have provided 
here a fourth perspective to look into GA and its local nature accompanying the infinite spatial 
dimension limit towards the GWF formalism. Furthermore, the proof readily suggests that such a 
simple scheme has superior performance than GA in finite systems. This fact is case studied here 
on small Hydrogen clusters. To distinguish it from GA, we call it A A in the coming discussions. 

COMPARISON BETWEEN THE NEW APPROXIMATION AND GA 

Let’s first look at the performance of AA and GA on predicting the ground state energy of 
a finite system constructed by Hydrogen atoms. We choose the one-dimensional minimum basis 
Hydrogen chains with periodic boundary conditions to carry out the calculation. The ab initio 
Hamiltonian is 

H = Uij^k,lclfjC^j^^iCk^cr'Cl,a (30) 

where (i,j) runs through all possible site pairs and {i,j,k,l) describes all possible 2-body inter¬ 
actions on the chain. The bare energy parameters t and U are evaluated with GAMESS (US), a 
Quantum Ghemistry package widely used for molecular calculations. Among all the interactions, 
the density-density interactions are the most dominant. Its ratio against nearest neighbor hopping 
is a quantitative measure of the strength of correlation, which increases as inter-atomic distance in¬ 
creases. The Gutzwiller wavefunction is constructed by applying the Gutzwiller projector defined 
in Eq. [3] onto the noninteracting Hartree-Fock wavefunction formed by the occupied molecular 
orbitals. We check the performance of both approximations on two systems formed by 8 and 16 
hydrogen atoms respectively. The brute force evaluation with the Gutzwiller trial wavefunction 
(GWF) and the Configuration Interaction (Cl) calculation give benchmark results for AA and GA 
to compare. As shown in Fig. [TJ GWF gives quite close energy to Cl and is able to reach the 
correct atomic energy given by CL One might note that Cl gives an atomic Hydrogen ground state 
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FIG. 1: Total energy evaluated using the current method (denoted as AA), GA, Coupled Cluster (CCSD), 
Configuration Interaction (Cl), and GWF if possible, on circular Hydrogen (H) chains of N=8 and N=16 
atoms. The ab initio Hamiltonians are returned by GAMESS(US) using a STO-3G basis set description of 
H. For N=8, GWF and Cl results are both available to benchmark AA and GA. For N=16, only CCSD 
is evaluated to check both approximations. Note CCSD does not converge well beyond the bond breaking 
region. GWF and Cl do not reach the known Hydrogen atomic energy of -O.SHartree due to the use of 
minimum basis set. The noninteracting wavefunction underlying GWF (thus GA and AA as well) are fixed 
at the Hartree-Fock solution. 


energy differring from the well-known value of —13.6eV. This is due to the STO-3G minimum basis 
set chosen to set up the ab initio Hamiltonian in Eq. [30j In the 16 Hydrogen atom system, GWF 
and Cl are not convenient to evaluate. Thus, we use the Coupled Cluster method(CCSD) instead 
to give an estimate on the ground state energy. CCSD is a very popular method in Quantum 
Chemistry ^ with good balance between speed and accuracy. It has some known issues but is 


good enough for our current purpose. The comparison between AA and GA, together with other 
methods, are presented in Fig. [U where the ground state energies are plotted against different 
inter-atomic separations. Several things are quite interesting to note. GA is in general not as 
good as AA, but its performance seems to be improved as system size increases. This gives a 
strong support on our proof of GA being a grand canonical ensemble description of the current 
coefficient averaging approximation, while AA is its canonical ensemble description. Another thing 
to notice is that, as compared with GWF, AA is slightly higher in energy around the bonding 
region and does not reach the atomic limit energy returned by GWF. This is mainly due to biased 
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FIG. 2: Shown in the figure are the renormalization factors defined in Eq. 1141 as a function of g for the 
onsite and nearest neighbor density-density interactions evaluated with GWF, AA and GA. The calculation 
is on a circular chain of 8 Hydrogen atoms described with a STO-3G basis set and separated with an atomic 
distance of 2.5A. The noninteracting wavefunction is provided by a restricted Hartree Fock calculation. The 
optimal Gutzwiller parameter has a value oi g ~ 0.02. Energy from AA is about 0.02 Har lower than GWF 
at r = 2.5A as seen from Fig[T] This energy difference is mainly due to the inaccurate (niqn 2 , 4 ,) evaluated 
with AA, the current approximation. 


onsite double occupancy as well as insufficient enhancement of off-site density-density correlations 
in the new approximation, as manifested in Fig [2] for an inter-atomic separation of 2.5A on the 
8 Hydrogen chain system. Consider that the optimal Gutzwiller parameter g is quite small, AA 
thus gives roughly consistent onsite energy as GWF, but underestimates the inter-site Coulomb 
repulsion. Consequently, AA gives an underestimated bounding energy. On the other hand, GA 
deviates more seriously from the benchmark energies. It also gives rise to an unphysical kink close 
to the bond breaking regime hinting a plausible meta-stable state for the Hydrogen rings. 

As the system size increases, renormalization factors, Zq, in the new approximation gradually 
approach the GA results, as illustrated in Fig. [3] and Fig. [4] on four typical operators in the 
half-filled case. Results deviating from half-filling are not shown as both methods are quite close to 
GWF results(e.g. see Ref for the nonhalf-filling behavior of GA). From Fig. [3] we immediately 
see noticeable finite size effects for AA. It converges quickly as system size increases. At a system 
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FIG. 3: Renormalization factor Z of four typical operators, c[^C 2 ,'\- and rii^an 2 ,a' are 

evaluated with AA for different linear H chain sizes. GA results are also provided for comparison. All graphs 
share the same legend given in graph (a). Note, no Hamiltonian is needed to scan g dependence of AA and 
GA, only the Gutzwiller wavefunction is relevant. 


size of = 10^, both AA and GA are nearly indistinguishable from each other. To provide a more 
quantitative measure on the asymptotic behavior of the new approximation, we present the system 
size dependence of the relative differences of renormalization factors between AA and GA in Fig. 
m The manifest linear dependence between these two quantities provides solid support on the 
conclusion that the new approximation asymptotically approaches GA in its infinite lattice limit. 
This might also suggest that the single band Gutzwiller approximation is quite an unexpectedly 
stable limit for the Gutzwiller wavefunction formalism with a site-wise local Gutzwiller projector. 

AN ATTEMPT TO APPLY THE NEW APPROXIMATION TO MULTIBAND SYSTEMS 

Thus far, we have applied the new approximation to single band systems and compared it 
against GA and other benchmark calculations. Encouraged by the limited success, we would like 
to generalize it to multiband systems which are physically more interesting and relevant to real 
problems. The simplest way to do so is to rewrite the Gutzwiller projector in terms of multi-orbital 
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FIG. 4: System size dependence of the relative renormalization factor differences between the new approx¬ 
imation and GA is shown for 3 different operators discussed in Fig. [31 The calculation is carried out with 
g = 0.2. C 2 j. “ C 21 both approximations and is thus suppressed from the plot. Note the log 

scale is used on both axes. 


particle number operators 


3.3l |. and to follow Eq. [T8]to Eq. [23]to impose overall physical constraints. 


Unfortunately, this is not successful due to the very fact that the proposed approximation itself does 
not automatically preserve the correct total number of electrons in the system except for the single 
band case. Indeed, there is no guarantee that relevant physical constraints can be automatically 
preserved with such a simple and arbitrary wavefunction coefficient averaging scheme. It is a bit 
of luck that it does work out for single band systems. 

However, this is not the excuse for us to give up the whole idea as it does provide us with a 
candidate recipe to define the renormalization factor for any operator. If we could look into the 
issue of failure a bit closer and at the same time fill us with confidence that the idea of averaging 
does work well on single band systems, we might come up with the following interpretation on the 
validity of the scheme: the averaging procedure does capture part of the correlation physics by giv¬ 
ing reduced hopping or onsite double occupancy, etc., but it might describe correlation physics only 
qualitatively right in generic cases. One way to improve its quality requires more degrees of freedom 
to enable adjusting the way how the averaging can be carried out. This is actually quite similar to 
the logic Bunemann et. al. followed in deriving the multiband Gutzwiller approximation^ 32]. 


There, a set of orbital-wise fugacity parameters are introduced into the Gutzwiller trial wavefunc¬ 
tion as adjustable parameters to ensure that the orbital occupations take the values evaluated with 
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the noninteracting wavefunction 


33j], or fulfill some constraints of similar nature in case the local 


Gutzwiller projector is more general than involving only particle operators 


3d | 38|. Their purpose 


to choose these constraints is for mathematical convenience to enable a rigorous Feynmann di¬ 
agrammatic expansion to the total energy expectation value evaluated with the Gutzwiller trial 
wavefunction. Similarly, we will introduce the fugacity parameters into the trial wavefunction and 
impose other more reasonable physical constraints than creating a connection with information 
from the noninteracting wavefunction. We will use these fugacity parameters to calculate the 
renormalization factors defined through averaging. From a statistical point of view, these fugacity 
parameters define a weight on each wavefunction coefficient. Thus the renormalization factors are 
defined not as a simple averaging, but as a weighted averaging instead. With a careful choice of 
weights, a concept readily borrowed from the statistics theory, more reasonable renormalization 
factors are anticipated to make this new approximation useful. Specifically, we introduce into the 
trial wavefunction of Eq. [2] a local weighting operator 


= (31) 

S 

where i is site index, s denotes each local correlated orbital and rjg is the weight of that state. Now 
the Gutzwiller trial wavefunction to be used in defining the renormalization factors is 



Note, this wavefunction is not supposed to be used to evaluate physical expectation values of 
any operator. They are evaluated with Eq. [2] directly, and are thus unique given a Gutzwiller 
parameter g. Of course, the choice of the weighting operator, Wi, is very flexible, depending on 
what kind of physical properties one decides to preserve. In case of a Hydrogen dimer described 
by a large basis set of 6-311G and to be studied here, we use Eq. [3T]to ensure correct total number 
of electrons in the dimer system. The correlated orbital is chosen to be Is while other orbitals 
are treated uncorrelated. With these enhanced definitions, the whole averaging formalism can be 
written out and r]s can be determined analytically as a function of g. Details of the calculation 
are given in Appendix B. Rather unexpectedly, such a slightly modified averaging scheme gives an 
exact description on the multiband Hydrogen dimer system within a Gutzwiller trial wavefunction. 
Besides the nice agreement, one might appreciate a piece of physics embedded in the current 
approximation. We found that the fugacity parameter has two solutions. One solution gives 
orbital occupations identical to those of the underlying noninteracting wavefunction, a scenario 
adopted in the multiband GA formalism, and the other gives the correct orbital occupation in 
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the Gutzwiller trial wavefunction. The coincidence of the solutions to the assumption underlying 
the multiband GA formalism might hint that the weighted version of the new approximation is 
intimately related with the multiband GA formalism that Bunemann et. al derived, just like AA to 
GA in the single band case. This fact might also speak loud that the new approximation introduced 
in this work could perform better than GA in capturing the strong correlation physics inherent in 
the Gutzwiller wavefunction formalism. 


CONCLUSION 


In sum, we have introduced a very simple yet effective approximation towards the Gutzwiller 
wavefunction formalism. The simplicity is easily seen by noting that the renormalization factor 
of an operator is obtained through a direct average over nonvanishing coefficients of a Gutzwiller 
wavefunction constructed with a density-density type Gutzwiller projector. The effectiveness is 
supported by the agreement between the new approximation and GA, a well-studied approxima¬ 
tion towards the Gutzwiller wavefunction formalism, in the infinite lattice limit of single band 
systems. Thus the current approximation provides a new perspective towards GA and its under¬ 
lying assumptions. The proof showing their mutual relationship in Appendix A clearly reveals the 
grand canonical ensemble nature of GA, which readily prevents it from being applied to finite sys- 


tems. This very nature of GA also shows itself in the counting argument interpretation of GA 


but is not clear from the diagram based formulation 


m 


33t |. In this work, several numerical instances 


are provided to compare the performance of the new approximation and GA, and to show the 
asymptotic agreement between these two schemes. 

Although it shows an improved performance than GA in single band systems, the naive aver¬ 
aging, however, does not capture the double occupancy and other subtle physical quantities well 
enough. It also leads to the failure to reach the correct local orbital occupation in a Hydrogen 
dimer described with a large basis set, the simplest multiband system to be studied. All these can 
be improved however, at least partially, by introducing more fugacity parameters and physical con¬ 
straints into the scheme. For example, a simple modification of introducing a weighting operator 
like Eq. [ST] nicely fixes the multi-orbital Hydrogen dimer problem. Actually, such a simple problem 
indicates a possible close relation between the weighted averaging scheme and the multiband GA 
formalism. With a more carefully designed weighting factor and more physical constraints included 
in the formalism, the performance of the current approximation could be systematically improved 
to reach that of the Gutzwiller wavefunction formalism, with which the current approximation is 
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built to match. 

This new idea is applied most conveniently towards Gutzwiller trial wavefunction with its 
Gutzwiller projector commuting with number operators. Such a density-density type Gutzwiller 
projector (d-GPJ) might sound inferior than a more generic Gutzwiller projector (g-GPJ) defined 


with a general interaction operator to describe local correlation effects [3^. This might cast doubt 
on the ultimate usefulness of the current scheme to describe real systems. We have to admit that 
more variational degrees of freedom introduced in g-GPJ are able to give ground state energy. But 
we also believe that d-GPJ might perform with close quality as g-GPJ. It might also have an addi¬ 
tional advantage of being more convenient in its implementation in practice. To fill us with some 
confidence first, we might find a successful application of d-GPJ in capturing the correct physics. 
Actually, such a Gutzwiller trial wavefunction has been applied to the Hubbard model and leads to 
the well-known t-J Hamiltonian in the strong correlation limit 29|]. The t-J Hamiltonian, as well as 


the Heisenberg Hamiltonian as its special case, is the main working horse towards strong correla¬ 
tion limit and spin dynamics. Second, we believe that g-GPJ can only be treated in the multiband 
GA level in practice. Thus a much fairer comparison would be between the multiband GA and the 
current scheme, to compare which one is able to perform better. Then, the answer is clear that 
no party wins over the other for sure. One motivation to develop the formalism involving g-GPJ 
might be because GA from d-GPJ has to assume {n(a,a )) = a) > ^ serious limitation in applying 
GA to multiband systems 331]. While GA from g-GPJ has this constraint removed, it still adopts 
constraints of similar nature in its derivation in order to carry out a Feynmann diagrammatic 
expansion 3^. Thus, how far multiband GA can be free from this constraint is unknown to us. In 


the new approximation introduced in this work, however, local orbital renormalization inherently 
exists which necessarily gives a local orbital occupation differring from its noninteracting value. 
One might systematically introduce more and more physical constraints into the formalism to help 
guide its outcome to be more and more physical. One more advantage of the current scheme against 
GA is that it is readily applied to Gutzwiller projectors with nonlocal density-density correlations. 
This fact is quite easy to notice although its practical implementation might be an issue. 
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tions. This work is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, 
Division of Materials Sciences and Engineering. Ames Laboratory is operated for the U.S. Depart¬ 
ment of Energy by Iowa State University under Gontract No. DE-AC02-07CHII358. 
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To establish a connection between these two approximations, let’s consider the renormalization 
coefficient for O {i,j) = cj^Cj^a- with site indices i / j. Here a denotes a spin-orbital composite 
state. Following the spirit leading to Eq. [18] to Eq. [20l one can write down and H for as 


Z^{r} Z]{r'} (fl'rifl'rjfl'r'fl'r'J S ^0;r,r'^ 

Z]{r} S{r'} ^ <^r 

Ein 

E(r) 4” 


(33) 


(34) 


Here {F} and {F'} denote a complete set of occupation conhguration on the whole lattice, 
denotes the g factor at site i with local Fock state Fj. In order to reach GA, let’s hrst relax the 8 
constraints, and note the following relations hold 



'^l = V^ 
{r} 


(35) 

(36) 


with N total number of lattice sites and T> dimensionality of the local Fock space. We thus have 


T 

^=B = 


(Srfi'r) 


E 


5ri5rj5r'5r' 


T,;,r 




(r.lTIU) |(r,|c,,qr')| 


(37) 


This expression for rj is only intermediate and not fully consistent, as one can see that rj fails to 
vanish as variational parameters g in the standard Gutzwiller projector approaches 0. 

Now let us introduce another approximation that the probability of hnding a given Fock space 
occupation is a product of probabilities of hnding the specihc occupation conhguration on each 
lattice site, or namely. 


l({r}|'l')l" = npr. (38) 

k 

|({Ffc}|To)|2 = nPr, (39) 

k 

where Pr^^Pr,. denote the probability to hnd a Fock state F^ on site k for the trial and nonin¬ 
teracting wavefunction respectively. This is too big a step forward as |({F} |'k)|^ is normally not 
decomposeable on lattice sites. This approximation can be validated only in the inhnite spatial 
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dimension limit 4^, which also underlies the Gutzwiller approximation. With this approximation, 
it is reasonable to assume 


STk = 


i PTk 

Pk 


(40) 


for the Gutzwiller trial wavefunction defined in Eq. [2] by noting that the local projection operator 
is site-wise. Feed Eq. 00] into Eq. [37] and we get 


ri = ZiZj 


with 


Zi = 


E 


PTiPr' 


{j:rPr/4) 


|(Ulc,,,|r')| 


(41) 


(42) 


A mean-field type approximation is now introduced for terms involving by replacing those 
terms with their averages. Again, in the infinite spatial dimension limit and under some convenience 
conditions, can be expressed in terms of local occupancy, on each local state a, as 


o-er cr^r 

Then there are two averages needed to be calculated. 


'' 2 / Ti.r' 

11 -I 


|r:>l 


TiY E 

V 2 ) o-^r' 


Ucr 


and 


= -\/n0 (1 -nO) 


r 


(43) 

(44) 

(45) 

(46) 

(47) 


Here T) is the dimensionality of the local Fock space. The prefactor 2 in Eq. 061 accounts for the 
fact that Tj must not contain state a in it. Plug Eq. 061 and Eq. 071 back to Eq. 021 and one 
recovers the standard GA definition of Zi given in Eq. [5j 


APPENDIX B: THE NEW APPROXIMATION WITH WEIGHTED AVERAGE APPLIED 

ON THE MULTIBAND H DIMER SYSTEM 


To make things simple, let’s consider a local basis set composed of one correlated orbital, denoted 
as s, and 94 uncorrelated orbitals to describe the H-dimer, all defined as Wannier functions such 
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that they are orthogonal to each other within and between sites. A molecular orbital is thus 
generically created via operator 

4 = hiscl^^ + ( 48 ) 

i 

where c|^ ^ creates an electron at site i and orbital a with spin a, or, c|^o-|0) = \ia,a). The 
coefficients satisfy the normalization condition 

Ihi.lV = l (49) 

OL^S 

with translational invariance implicitly assumed. That the cross terms contributing to the normal¬ 
ization vanish comes from the fact that each atomic orbital is a Wannier function, as mentioned at 
the beginning of this appendix. Eq. 09] also expresses the electron conservation condition, ensuring 
each site has half an electron with a specific spin. The noninteracting Hartree-Fock wavefunction 
can be expressed as 

|To) = a|a||0) (50) 

for the H dimer, and the Gutzwiller wavefunction is defined as 
|\J/^ = I^O) 

= {g-l)\hu\^ |ls t, Is i) + (5 - 1) \h 2 s? |2s t, 2s i) + |To) (51) 


with g the Gutzwiller parameter. Expectation value w.r.t the Gutzwiller wavefunction for any 
operator can be straightforwardly evaluated. Specifically, for particle occupations, there are 


l)\hlsf + l I |2 

^ 2{g^-l)\huf + l' 

(^lQ!,cr) = - „ . I , ,4 l^lal for a ^ s 

2{g^-l)\hu\ -hi 


(52) 

(53) 


Obviously, it is easy to test that they satisfy the electron conservation condition, Eq. 09] with help 
of that constraint. 

For the current approximation with weighted average enhancement, the local weighting operator 
is chosen as 
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The expectation value for an operator is considered from its denominator, wavefunction normal¬ 
ization, and its numerator respectively. For the normalization factor (T|T), there is 


^ + 2r?^ + + (2C|g^ + 2f)T) 

2774 + 2ry4 + 4(^2^^ + (2C2^ + 2 fn) ^ 

with C™ the usual combinatorial number choosing m elements out of n elements. In the numer¬ 
ator of the renormalization factor in Eq. IMl each term has clear physical interpretation, 
corresponds to two electrons occupying s orbitals on the same site, has the two electrons take s 
orbitals on different sites, is related to Fock states with only one electron in s orbitals, while the 
last term corresponds to Fock states with no s orbital. The terms in the denominator of the above 
prefactor are obtained by ignoring the variational parameter g. What is left acts as the weight to 
each term in the numerator, a necessary step to normalize a weighted average. Similarly, one can 
write down the numerator of the expectation value of an operator. For electron occupations, they 
are 



(^1 

(^1 ^la,t 


|T) ~ 


q^g"^ + q^ + 

2q* -h 


|T) ~ \hia\‘^ for a 7 ^ s 


I his 


(55) 

(56) 


The constraint of conserved electron occupation on a H dimer requires 

(T| his,^ |T) + ^ (T| hia,t 1^) = \ (^1^) (57) 

Feed Eq. [SH Eq. [55] and Eq. [56] into Eq. [57] and solve for q, and one ends up with two solutions 


77 ^ = 0 


77 ^ = 


011/7 


Is 


(58) 

(59) 


\hla\ 

Interestingly, both these two solutions have clear physical interpretations. 77 = 0 corresponds to 
the case where the correlated system has the same local orbital occupation as the underlying non¬ 
interacting wavefunction, which is what Bunemann’s multiband Gutzwiller approximation starts 
with. The nontrivial solution of Eq. [59] to 77 gives the correct charge occupation as the rigorous 
Gutzwiller wavefunction. One can readily verify this fact by inserting the solution Eq. [59] back to 
the expressions for ( 77 is,'f) and and take the constraint of Eq. [l9] Actually, one can further 

verify that this nontrivial solution renders correct expressions for any one and two body operators 
of the H dimer system. 
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